function post_stress_T
fold=fopen('153.bdf','r');
fnew=fopen('153+.bdf','w');
while feof(fold)==0
    entry=fgetl(fold);
    if isempty(entry)
        fprintf(fnew,'\n');
        continue
    end
    for i=1:min(5,length(entry))
        if entry(i)=='$'
            entry(i)=[];
            fprintf(fnew,'%s\n',entry);
            fgetl(fold)
            fgetl(fold)
            entry=fgetl(fold);
        end
    end
    fprintf(fnew,'%s\n',entry);
end
fclose('all');
!del 153.bdf
!rename 153+.bdf 153.bdf
if ~exist('plot','dir')
    mkdir('plot');
end
load result\T_up_596_600.mat T_up
load result\T_down_596_600.mat T_down
load temp\str_data.mat str_data
load temp\T.mat T
load temp\W.mat W
load temp\stru_grid.mat stru_grid
load temp\P_aero_up.mat P_aero_up
load temp\P_aero_down.mat P_aero_down
load temp\Q_up.mat Q_up
load temp\Q_down.mat Q_down
load temp\D_T.mat D_T
stru_grid_out=stru_grid((stru_grid(:,1)>400000)|((200000<stru_grid(:,1))&(stru_grid(:,1)<300000)),:);
W=W(:,end); W=reshape(W,6,size(W,1)/6)';
stru_grid=str_data.stru_grid; aerobox=str_data.aerobox;
stru_up_rigid=str_data.stru_up_rigid; stru_down_rigid=str_data.stru_down_rigid;
stru_quad_up=str_data.stru_quad_up; stru_quad_down=str_data.stru_quad_down;
stru_quad_up_plus=str_data.stru_quad_up_plus; stru_quad_down_plus=str_data.stru_quad_down_plus;
[aero_up_pre,aero_down_pre,~,~,~,~]=interpolate_point(1000,aerobox,stru_quad_up);
Npath='E:\nastran\bin\nast20180.exe';
dt=4;
t_end=600;
try
    load temp\stress_T.mat stress_T
catch
    stress_T=[];
end

for ii=size(stress_T,2)*4:dt:t_end
    try
        copyfile(['result\159_',num2str(ii),'_',num2str(ii+dt),'.pch'],'159.pch');
        fold=fopen('159.pch','r');
        fnew=fopen('159+.pch','w');
        while feof(fold)==0
            entry=fgetl(fold);
            if length(entry)<2
                continue
            end
            if length(strfind(entry,'TEMP*                  1'))==1
                entry(44:56)=' 3.000000E+02';
            end
            fprintf(fnew,'%s\n',entry);
        end
        fclose('all');
        !del 159.pch
        !rename 159+.pch 159.pch
        run_153(Npath)
        !del 153.f04;
        !del 153.f06;
        !del 153.log;
        !del 153.MASTER;
        !del 153.DBALL;
        !del 153.IFPDAT;
    
        f=fopen('153.pch');
        z=1;
        stress=zeros(0,2);
        while true
            tline1=fgetl(f);
            if length(strfind(tline1,'$LOAD FACTOR =   1.0000000E+00'))==1
                while true
                    tline1=fgetl(f);
                    if tline1(1)=='$'
                        break
                    end
                    stress(z,1)=sscanf(tline1(1:16),'%d');
                    for j=1:3
                        tline1=fgetl(f);
                    end
                    stress(z,2)=sscanf(tline1(41:58),'%f');
                    z=z+1;
                    for k=1:45
                        fgetl(f);
                    end
                end
                break
            end
        end
        stress(:,2)=stress(:,2)/1000;
        f=fopen(['plot\stress_T',num2str(ii),'_',num2str(ii+dt),'.plt'],'w');
        fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
            ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
            ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
        l=1;
        for i=1:size(stru_grid,1)
            fprintf(f,'%e ',stru_grid(i,2));
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        fprintf(f,'\n');
    
        l=1;
        for i=1:size(stru_grid,1)
            fprintf(f,'%e ',stru_grid(i,3));
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        fprintf(f,'\n');
    
        l=1;
        for i=1:size(stru_grid,1)
            fprintf(f,'%e ',stru_grid(i,4));
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        fprintf(f,'\n');
    
        l=1;
        for i=1:size(stress,1)
            if ismember(stress(i,1)-200000,stru_quad_up)
                fprintf(f,'%e ',stress(i,2));
            else
                continue
            end
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        for i=1:size(stress,1)
            if ismember(stress(i,1)-200000,stru_quad_down)
                fprintf(f,'%e ',stress(i,2));
            else
                continue
            end
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        for i=1:size(stress,1)
            if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
                fprintf(f,'%e ',stress(i,2));
            else
                continue
            end
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        for i=1:size(stress,1)
            if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
                fprintf(f,'%e ',stress(i,2));
            else
                continue
            end
            if l==5
                fprintf(f,'\n');
                l=1;
            else
                l=l+1;
            end
        end
        fprintf(f,'\n');
    
        p=0;
        for i=1:length(P_aero_up)
            fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
                find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
                find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
                find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
            p=p+1;
        end
        for i=1:length(P_aero_down)
            fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
                find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
                find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
                find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
            p=p+1;
        end
        for i=1:size(stru_quad_up_plus,1)
            if stru_quad_up_plus(i,1)==-1
                break
            end
            fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
                find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
                find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
                find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
            p=p+1;
        end
        for i=1:size(stru_quad_down_plus,1)
            if stru_quad_down_plus(i,1)==-1
                break
            end
            fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
                find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
                find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
                find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
            p=p+1;
        end
        fclose('all');
        !del 159.pch;
        !del 153.pch;
        stress_T(:,ii/4+1)=stress(:,2);
    catch
        stress_T(:,ii/4+1)=0;
        files_ini;
    end
    save temp\stress_T.mat stress_T
end
stress_T_max=zeros(151,2);
stress_T_max(:,1)=1:151;
for i=1:size(stress_T,1)
    stress_T_max(i,2)=max(stress_T(:,i));
end
save temp\stress_T_max.mat stress_T_max

stress_T_max(stress_T_max(:)==0,:)=[];
stress_T_max(stress_T_max(:)>450,:)=[];
plot(stress_T_max(:,1),stress_T_max(:,2))
end